#*************        丁香园-NHANES     *************
#*************        Analyze Code      *************
#*************        倾向评分讲解      *************
#*************                          *************
# 倾向评分匹配是一种统计学方法，用于处理观察研究的数据。
# 在观察性研究中（NHANES 是典型的观察性研究），由于种种原因，数据偏差和混杂变量较多。
# 倾向评分匹配的方法正是为了减少这些偏差和混杂变量的影响，
# 以便对核心研究变量 X 的不同取值组进行更合理的比较。

#### 0.准备好环境 ####
library(gtsummary)
library(survey)
library(haven)
library(tableone)
library(rms) # RCS-Regression Modeling Strategies
library(MatchIt) # PSM
library(plyr)
library(tidyverse)
setwd('G:/BaiduNetdiskDownload/NHANES')
# 运行一下读取 mort 数据的函数
NHSMortRead <- function(file.path){
  srvyin <- file.path   # full .DAT name here
  srvyout <- "data" # shorthand dataset name here
  
  # Example syntax:
  #srvyin <- paste("NHANES_1999_2000_MORT_2019_PUBLIC.dat")   
  #srvyout <- "NHANES_1999_2000"      
  
  
  # read in the fixed-width format ASCII file
  dsn <- read_fwf(file=srvyin,
                  col_types = "iiiiiiii",
                  fwf_cols(seqn = c(1,6),
                           eligstat = c(15,15),
                           mortstat = c(16,16),
                           ucod_leading = c(17,19),
                           diabetes = c(20,20),
                           hyperten = c(21,21),
                           permth_int = c(43,45),
                           permth_exm = c(46,48)
                  ),
                  na = c("", ".")
  )
  
  # NOTE:   SEQN is the unique ID for NHANES.
  
  # Structure and contents of data
  str(dsn)
  
  
  # Variable frequencies
  
  #ELIGSTAT: Eligibility Status for Mortality Follow-up
  table(dsn$eligstat)
  #1 = "Eligible"
  #2 = "Under age 18, not available for public release"
  #3 = "Ineligible"
  
  #MORTSTAT: Final Mortality Status
  table(dsn$mortstat, useNA="ifany")
  # 0 = Assumed alive
  # 1 = Assumed deceased
  # <NA> = Ineligible or under age 18
  
  #UCOD_LEADING: Underlying Cause of Death: Recode
  table(dsn$ucod_leading, useNA="ifany")
  # 1 = Diseases of heart (I00-I09, I11, I13, I20-I51)
  # 2 = Malignant neoplasms (C00-C97)
  # 3 = Chronic lower respiratory diseases (J40-J47)
  # 4 = Accidents (unintentional injuries) (V01-X59, Y85-Y86)
  # 5 = Cerebrovascular diseases (I60-I69)
  # 6 = Alzheimer's disease (G30)
  # 7 = Diabetes mellitus (E10-E14)
  # 8 = Influenza and pneumonia (J09-J18)
  # 9 = Nephritis, nephrotic syndrome and nephrosis (N00-N07, N17-N19, N25-N27)
  # 10 = All other causes (residual)
  # <NA> = Ineligible, under age 18, assumed alive, or no cause of death data available
  
  #DIABETES: Diabetes Flag from Multiple Cause of Death (MCOD)
  table(dsn$diabetes, useNA="ifany")
  # 0 = No - Condition not listed as a multiple cause of death
  # 1 = Yes - Condition listed as a multiple cause of death
  # <NA> = Assumed alive, under age 18, ineligible for mortality follow-up, or MCOD not available
  
  #HYPERTEN: Hypertension Flag from Multiple Cause of Death (MCOD)
  table(dsn$hyperten, useNA="ifany")
  # 0 = No - Condition not listed as a multiple cause of death
  # 1 = Yes - Condition listed as a multiple cause of death
  # <NA> = Assumed alive, under age 18, ineligible for mortality follow-up, or MCOD not available
  
  # Re-name the dataset, DSN, to the short survey name then remove other R objects
  assign(paste0(srvyout), dsn)
  rm(dsn, srvyin, srvyout)
  # View(data)
  return(data)
}


#### 1. 倾向评分（加权情况下） ####
##### 1.1 提取数据模块 ##### 
demo.i <- read_xpt("2015-2016/Demographics/demo_i.xpt")#要提前设置好数据存储的路径
colnames(demo.i)
dim(demo.i) # 9971

# 提取生存数据 ：https://www.cdc.gov/nchs/data-linkage/mortality-public.htm

# 连接 National Death Index (NDI) ，补充 NHANES 参与人员的死亡数据
# 得到公开的 linked mortality files (LMF) 提供从调查参与日期到2019年12月31日的死亡率随访数据。
# 公开的数据中仅有随访时间或潜在的死亡原因
# 生存地址原始下载链接：https://ftp.cdc.gov/pub/Health_Statistics/NCHS/datalinkage/linked_mortality/
# 可以下载不同周期的生存数据，下载后进行读取

# setwd("C:/PUBLIC USE DATA")
mort.data_2015_2016 <- NHSMortRead('NHANES_2015_2016_MORT_2019_PUBLIC.dat')
# eligstat
# A value of 1 for ELIGSTAT indicates that the survey participant was eligible for the mortality linkage, 
# a value of 2 indicates the survey participant was under age 18 and not eligible for public release, 
# a value of 3 indicates the survey participant was not linkage eligible due to having insufficient identifying data to conduct data linkage.

# ucod_leading: 10种 潜在死亡原因，001-010，分别对应的原因看课程附件材料
# mortstat: 1生存状态
# permth_int: 从interview到随访时间点，间隔的月份数
# permth_exm: 从MEC-移动检查车到随访时间点，间隔的月份数

dim(mort.data_2015_2016) # 9971
# 转变为大写的列名，和 demo 进行拼接
colnames(mort.data_2015_2016) <- toupper(colnames(mort.data_2015_2016))
# View(mort.data_2015_2016)


##### 1.2 合并数据提取变量  & 去掉 NA 的行 & 衍生变量 ##### 
analyze.sample.data <- demo.i[,c('SEQN',
                                 "RIDAGEYR", "RIAGENDR", "INDFMPIR", "DMDEDUC2",
                                 "WTINT2YR", "SDMVPSU", "SDMVSTRA")]

analyze.mort.data <- mort.data_2015_2016[,c('SEQN', 'ELIGSTAT', 
                                            "MORTSTAT", "PERMTH_INT")]
# 拼接 demo & mort 数据
analyze.sample.data.add.mort <- join_all(list(analyze.sample.data, analyze.mort.data))

dim(analyze.sample.data.add.mort) # 9971    
# View(analyze.sample.data.add.mort) 

# 一键去掉 NA 的行
analyze.sample.data.add.mort.drop.na <- drop_na(analyze.sample.data.add.mort) 
dim(analyze.sample.data.add.mort.drop.na) # 5071

# 衍生低收入 vs 非低收入的变量 PIR.factor
analyze.sample.data.add.mort.drop.na$PIR.factor <- ifelse(analyze.sample.data.add.mort.drop.na$INDFMPIR < 1.3, 1, 0)

# 转换为分类变量
analyze.sample.data.add.mort.drop.na$DMDEDUC2 <- factor(analyze.sample.data.add.mort.drop.na$DMDEDUC2)

##### 1.3 生成复杂抽样 NHANES_design ##### 
NHANES_design <- svydesign(
  data = analyze.sample.data.add.mort.drop.na, 
  ids = ~SDMVPSU, 
  strata = ~SDMVSTRA, 
  nest = TRUE, 
  weights = ~WTINT2YR,
  survey.lonely.psu = "adjust") # 可以加上 survey.lonely.psu = "adjust" 避免1个PSU报错

summary(NHANES_design)

##### 1.4 判断是否需要进行 PSM：Table1 表格，不同 PIR.factor 下性别、年龄、教育成都的差异是否显著 #####
# 复杂抽样下的 Table1
tbl_svysummary(NHANES_design, by = PIR.factor, missing = 'no',
               include = c(RIDAGEYR, RIAGENDR, DMDEDUC2))%>%    
  add_p() 

# 替换 N 为非加权情况
tbl_svysummary(NHANES_design, by = PIR.factor, missing = 'no',
               include = c(RIDAGEYR, RIAGENDR, DMDEDUC2), # Use include to select variables
               statistic = list(all_continuous()  ~ "{mean} ({sd})",
                                # all_categorical() ~ "{n}    ({p}%)"),
                                all_categorical() ~ "{n_unweighted} ({p}%)"))%>% # 替换为非加权的 n
  modify_header(all_stat_cols() ~ "**{level}**, N = {n_unweighted} ({style_percent(p)}%)")%>%    
  add_p() 
  
  
##### 1.5 设置复杂抽样下的 PSM 权重 #####
ori.weight <- 1/(NHANES_design$prob)
mean.weight <- mean(ori.weight)
new.weights <- ori.weight/mean.weight # 生成用于计算的权重变量，加入到 data 中
data.for.machit <- cbind(NHANES_design$variables, new.weights) 
dim(data.for.machit) # 5071   13
# View(data.for.machit)

# 使用默认参数

 design.machit <- matchit(PIR.factor ~ RIDAGEYR + RIAGENDR + factor(DMDEDUC2), 
                          data = data.for.machit, s.weights = ~ new.weights)
# 使用较好的参数组合
design.machit <- matchit(PIR.factor ~ RIDAGEYR + RIAGENDR + factor(DMDEDUC2), 
                         data = data.for.machit, s.weights = ~ new.weights,
                         method = "nearest",
                         distance = "logit",
                         replace = FALSE,
                         caliper = 0.05)
# 得分差异上限（caliper）：当我们匹配用户的时候，我们要求每一对用户的得分差异不超过指定
# 的 caliper，“强扭的瓜不甜”，匹配不好就不要放弃吧。by 某一网友

# 查看函数的核心参数
# help(matchit)

# 查看 Standarized Mean Difference (SMD)
# 一般如果一个变量的 SMD 不超过 0.1，一般就可以认为这个变量的配平质量可以接受。
# 当一个变量的 SMD 超过 0.1 的时候，需要凭经验确认一下那个变量是不是没有那么重要。
# 我们还可以用倾向性得分来对用户进行加权，称为 Inverse Propensity Score Weighting (IPSW)。

summary(design.machit, standardize = T)$sum.matched

# 得到匹配后的数据
match.data.new <- match.data(design.machit)
# View(match.data.new)
dim(match.data.new) # 3282   16

##### 1.6 更新 design #####  
NHANES_design.after.psm <- svydesign(
  data = match.data.new, 
  ids = ~SDMVPSU, 
  strata = ~SDMVSTRA, 
  nest = TRUE, 
  weights = ~new.weights, # 更新 weight 的结果
  survey.lonely.psu = "adjust") # 可以加上 survey.lonely.psu = "adjust" 避免1个PSU报错

summary(NHANES_design)

# 得到 PSM 校正后的结果
tbl_svysummary(NHANES_design, by = PIR.factor, missing = 'no',percent = "row", 
               include = c(RIDAGEYR, RIAGENDR, DMDEDUC2))%>%    
  add_p() 


tbl_svysummary(NHANES_design.after.psm, by = PIR.factor, missing = 'no',percent = "row", 
               include = c(RIDAGEYR, RIAGENDR, DMDEDUC2))%>%    
  add_p() 

##### 1.7 复杂抽样的 Cox 回归模型 ##### 
# 死亡和贫困指数之间的相关性
# 校正前
fit.before.psm <- svycoxph(Surv(PERMTH_INT, MORTSTAT) ~ PIR.factor + RIDAGEYR + RIAGENDR + DMDEDUC2, 
                           design = NHANES_design)
fit.before.psm
# 校正后
fit.after.psm <- svycoxph(Surv(PERMTH_INT, MORTSTAT) ~ PIR.factor + RIDAGEYR + RIAGENDR + DMDEDUC2, 
                          design = NHANES_design.after.psm)
# 注意年龄和性别的 P value 的变化
fit.after.psm

